﻿#include "MeshTool.h"
#include "Algorithm.h"
#include "Mesh.h"
#include "TopoTables.h"
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/ConvertNumComponentsToOffsets.h>
#include <vtkm/cont/DataSetBuilderUniform.h>
//#include <vtkm/worklet/Triangulate.h>
#include <vtkm/filter/geometry_refinement/worklet/Triangulate.h>
#include <vtkm/worklet/WorkletCellNeighborhood.h>

namespace MeshTool
{
namespace detail
{
struct NumFacesPerCell : WorkletVisitCellsWithPoints
{
  using ControlSignature = void(CellSetIn cell_set, FieldOut numFacesInCell);
  using ExecutionSignature = void(CellShape, _2);

  template<typename CellShapeTag>
  VTKM_EXEC void operator()(CellShapeTag shape, IdComponent& num_faces) const
  {
    CellFaceNumberOfFaces(shape, num_faces);
  }
};

struct BuildDuplicateFace : WorkletVisitCellsWithPoints
{
  using ControlSignature = void(CellSetIn cellset, FieldOut face_id, FieldOut face_info);
  using ExecutionSignature = void(CellShape, PointIndices, InputIndex, VisitIndex, _2, _3);
  using ScatterType = ScatterCounting;

  template<typename CellShapeTag, typename C2NType, typename FaceInfo>
  VTKM_EXEC void operator()(CellShapeTag cell_shape,
                            const C2NType& c2n,
                            Id inputIndex,
                            IdComponent visitIndex,
                            Id3& faceId,
                            FaceInfo& face) const
  {
    CellFaceCanonicalId(visitIndex, cell_shape, c2n, faceId);
    face = { inputIndex, static_cast<UInt8>(visitIndex) };
  }
};

struct BuildDuplicateFace1 : WorkletVisitCellsWithPoints
{
  using ControlSignature = void(CellSetIn cellset,
                                FieldOut face_id,
                                FieldOut face_info,
                                FieldOut num_faces);
  using ExecutionSignature = void(CellShape, PointIndices, InputIndex, _2, _3, _4);
  //using ScatterType = ScatterCounting;

  template<typename CellShapeTag, typename C2NType, typename FaceInfo>
  VTKM_EXEC void operator()(CellShapeTag cell_shape,
                            const C2NType& c2n,
                            Id cell_id,
                            Id3& faceId,
                            FaceInfo& face,
                            IdComponent& num_faces) const
  {
    CellFaceNumberOfFaces(cell_shape, num_faces);
    for (IdComponent side = 0; side < num_faces; side++)
    {
      CellFaceCanonicalId(side, cell_shape, c2n, faceId);
      face = { cell_id, static_cast<UInt8>(side) };
    }
  }
};

struct BuildUniqueFace : public WorkletReduceByKey
{
  using ControlSignature = void(KeysIn faceKey,
                                ValuesIn faceInfo,
                                //ReducedValuesOut face_link,
                                WholeArrayOut cell_links);
  using ExecutionSignature = void(_2, /* _3,*/ _3);

  template<typename FaceInfoPair, typename CellToCellType /*, typename FaceInfo*/>
  VTKM_EXEC void operator()(const FaceInfoPair& face_info,
                            //Vec<FaceInfo, 2>& face_link,
                            CellToCellType& cell_link) const
  {
    if (face_info.GetNumberOfComponents() == 2) // 内部面
    {
      //face_link[0] = face_info[0];
      //face_link[1] = face_info[1];

      auto cell_id = face_info[0].Get().first;
      auto face_id = face_info[0].Get().second;
      cell_link.Get(cell_id)[face_id] = face_info[1];

      cell_id = face_info[1].Get().first;
      face_id = face_info[1].Get().second;
      cell_link.Get(cell_id)[face_id] = face_info[0];
    }
    else if (face_info.GetNumberOfComponents() == 1) // 边界面
    {
      //face_link[0] = face_info[0];
      // TODO: 引入边界条件
      //face_link[1] = { -1, 0 };

      auto cell_id = face_info[0].Get().first;
      auto face_id = face_info[0].Get().second;
      cell_link.Get(cell_id)[face_id] = { -1, 0 };
    }
    else
    {
      this->RaiseError("face_duplicate.GetNumberOfComponents error.");
    }
  }
};

struct ComputeTopolgy : WorkletVisitCellsWithPoints
{
  using ControlSignature = void(CellSetIn, FieldInPoint, FieldOut, FieldOut);
  using ExecutionSignature = void(_2, _3, _4);

  template<typename PType>
  VTKM_EXEC void operator()(const PType& pts, vtkm::Vec3f& centroid, vtkm::Vec3f& spacing) const
  {
    if (pts.GetNumberOfComponents() == 4)
    {
      centroid = (pts[0] + pts[1] + pts[2] + pts[3]) / 4;
      spacing[0] = (pts[1] - pts[0])[0];
      spacing[1] = (pts[2] - pts[1])[1];
      spacing[2] = 1;
    }
    else if (pts.GetNumberOfComponents() == 8)
    {
      centroid = (pts[0] + pts[1] + pts[2] + pts[3] + pts[4] + pts[5] + pts[6] + pts[7]) / 8;
      spacing[0] = (pts[1] - pts[0])[0];
      spacing[1] = (pts[2] - pts[1])[1];
      spacing[2] = (pts[4] - pts[0])[2];
    }
    else
    {
      RaiseError("unfinished");
    }
  }
};

struct BuildBCMape : vtkm::worklet::WorkletCellNeighborhood
{
  using ControlSignature = void(CellSetIn, FieldOut bnd_id);
  using ExecutionSignature = void(Boundary, _2);

  template<typename BndIdType>
  VTKM_EXEC void operator()(const vtkm::exec::BoundaryState& boundary, BndIdType& bnd_id) const
  {
    bnd_id = { -1, -1, -1, -1 };
    if (boundary.IJK[0] == 0) // left, inlet
    {
      bnd_id[3] = 3;
    }
    if (boundary.IJK[0] == boundary.PointDimensions[0] - 1) // right, outlet
    {
      bnd_id[1] = 1;
    }
    if (boundary.IJK[1] == boundary.PointDimensions[1] - 1) // top, wall
    {
      bnd_id[2] = 2;
    }
    if (boundary.IJK[1] == 0) // bottom
    {
      bnd_id[0] = 0;
    }
  }
};

struct BuildBCType : vtkm::worklet::WorkletMapField
{
  using ControlSignature = void(FieldIn element_side, FieldIn bc, WholeArrayOut cell_links);
  using ExecutionSignature = void(_1, _2, _3);

  template<typename FaceInfoPair, typename CellToCellType>
  VTKM_EXEC void operator()(const FaceInfoPair& face_info,
                            UInt8 bc_type,
                            CellToCellType& cell_link) const
  {
    auto cell_id = face_info.first;
    auto face_id = face_info.second;
    auto& element_side = cell_link.Get(cell_id);
    //std::cout << int(face_id) << std::endl;
    //std::cout << "N: " << element_side.GetNumberOfComponents() << std::endl;
    element_side[face_id] = { -1, bc_type };
    //std::cout << element_side[face_id].Get().first << std::endl;
    //[face_id];
    //std::cout << element_side[face_id].Get().first << std::endl;
    //std::cout << /*static_cast<int>*/(element_side[face_id].Get().second) << std::endl;

    //std::cout << cell_id << std::endl;
    //std::cout << face_id << std::endl;
  }
};

}

struct UpdateCoordinateSyatemWorklet : vtkm::worklet::WorkletMapField
{
  using ControlSignature = void(FieldIn in, FieldOut out);
  using ExecutionSignature = void(_1, _2);
  template<typename InType, typename OutType>
  VTKM_EXEC void operator()(const InType& in, OutType& out) const
  {
    out = in;
  }
};


void BuildFaceTopo(Mesh& mesh, const std::map<std::string, UInt8>& bnd_name_id)
{
  using FaceInfo = typename Mesh::FaceInfo;
  ArrayHandle<IdComponent> num_face_per_cell;
  Invoker{}(detail::NumFacesPerCell{}, mesh.GetCellSet(), num_face_per_cell);

  ArrayHandle<Id3> face_id;
  ArrayHandle<FaceInfo> face_info; //包含重复face
  //Invoker{}(
  //  detail::BuildDuplicateFace1{}, mesh.GetCellSet(), face_id, face_info, num_face_per_cell);

  Invoker{}(detail::BuildDuplicateFace{},
            ScatterCounting{ num_face_per_cell },
            mesh.GetCellSet(),
            face_id,
            face_info);

  auto& face_num_offset = mesh.GetFaceNumOffset();
  mesh.GetCellLinks().Allocate(face_info.GetNumberOfValues());
  face_num_offset = vtkm::cont::ConvertNumComponentsToOffsets(num_face_per_cell);
  num_face_per_cell.ReleaseResourcesExecution();

  auto cell_links =
    vtkm::cont::make_ArrayHandleGroupVecVariable(mesh.GetCellLinks(), face_num_offset);
  Keys<Id3> face_keys(face_id);
  face_id.ReleaseResourcesExecution();
  Invoker{}(detail::BuildUniqueFace{}, face_keys, face_info, cell_links);
  face_info.ReleaseResourcesExecution();




  //ArrayHandle<Vec<vtkm::Int16, 4>> bnd_id;
  //Invoker{}(detail::BuildBCMape{}, mesh.GetCellSet(), bnd_id);
  //auto& portal = bnd_id.ReadPortal();

  //for (size_t element = 0; element < bnd_id.GetNumberOfValues(); element++)
  //{
  //  auto& bc_type = portal.Get(element);
  //  for (size_t side = 0; side < 4; side++)
  //  {
  //    auto bc_id = bc_type[side];
  //    if (bc_id >= 0)
  //    {
  //      element_side.emplace_back(element, side);
  //      bc.push_back(bc_id);
  //    }
  //  }
  //}

  std::vector<FaceInfo> element_side;
  std::vector<UInt8> bc_type;

  auto valid_sideset_name = [](decltype(mesh.GetElementSide()) bc_element_side) -> std::string
  {
    std::stringstream ss;
    ss << "[";
    for (auto& var : bc_element_side)
    {
      ss << var.first << ",";
    }
    ss << "]";
    return ss.str();
  };

  auto& bc_element_side = mesh.GetElementSide();

  for (auto& var : bnd_name_id)
  {
    auto& sideset_name = var.first;
    try
    {
      void(bc_element_side.at(sideset_name));
    }
    catch (const std::exception&)
    {
      Error("side set:",
            sideset_name,
            ",在网格不存在,可用的side sets: ",
            valid_sideset_name(bc_element_side));
    }

    auto& ele_side = bc_element_side[sideset_name];
    for (auto& face_info : ele_side)
    {
      element_side.push_back(face_info);
      bc_type.push_back(bnd_name_id.at(sideset_name));
    }
  }

  //auto bounary_info = mesh.FlatBoundaryInfo(bnd_name_id);
  //ArrayHandle<UInt8> bc_array =bounary_info.first ;
  //ArrayHandle<FaceInfo> element_side_array = bounary_info.second;
  ArrayHandle<FaceInfo> element_side_array = vtkm::cont::make_ArrayHandle(element_side);
  ArrayHandle<UInt8> bc_array = vtkm::cont::make_ArrayHandle(bc_type);

  std::cout << bc_array.GetNumberOfValues() << std::endl;
  std::cout << element_side_array.GetNumberOfValues() << std::endl;
  Invoker{}(detail::BuildBCType{}, element_side_array, bc_array, cell_links);
}

vtkm::cont::DataSet CreateStructuredMesh(const Id2& point_dim,
                                         const Vec2f& orign,
                                         const Vec2f& space)
{
  vtkm::cont::DataSet data_set = vtkm::cont::DataSetBuilderUniform::Create(point_dim, orign, space);

  ArrayHandle<Vec3f> spacing, center;
  //ArrayHandle<Id> b_id;
  //ArrayHandle<vtkm::Int16> blk_id;

  vtkm::cont::Invoker{}(detail::ComputeTopolgy{},
                        data_set.GetCellSet(),
                        data_set.GetCoordinateSystem(),
                        center,
                        spacing);

  data_set.AddCellField("spacing", spacing);
  data_set.AddCellField("center", center);
  //data_set.AddPointField("boundary_id", b_id);
  //data_set.AddCellField("blk_id", blk_id);

  return data_set;
}

vtkm::cont::DataSet CreateStructuredMesh(const Id3& point_dim,
                                         const Vec3f& orign,
                                         const Vec3f& space)
{
  vtkm::cont::DataSet data_set = vtkm::cont::DataSetBuilderUniform::Create(point_dim, orign, space);

  //ArrayHandle<Vec3f> spacing, center;
  //ArrayHandle<Id> b_id;
  ////ArrayHandle<vtkm::Int16> blk_id;

  //vtkm::cont::Invoker{}(detail::ComputeTopolgy{},
  //                      data_set.GetCellSet(),
  //                      data_set.GetCoordinateSystem(),
  //                      center,
  //                      spacing);


  //data_set.AddCellField("spacing", spacing);
  //data_set.AddCellField("center", center);

  return data_set;
}

void UpdateCoordinateSyatem(vtkm::cont::DataSet& data_set) 
{
  ArrayHandle<Vec3f> pts;
  Invoker{}(
    UpdateCoordinateSyatemWorklet{}, data_set.GetCoordinateSystem().GetDataAsMultiplexer(), pts);
  data_set.GetCoordinateSystem() = CoordinateSystem("coord", pts);
}


void TriangulateMesh(Mesh& mesh)
{
  vtkm::worklet::Triangulate triangulate;

  auto out_cellset = triangulate.Run(mesh.GetCellSet());
  mesh.GetCellSet() = out_cellset;
  //this->_cell_set = output.GetCellSet();
  //ArrayHandle<Vec3f> pts;
  //Invoker{}(CoodConvert{}, output.GetCoordinateSystem(), pts);
  //this->_coord = CoordinateSystem("coord", pts);
}

}